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This paper introduces a new structural phase field crystal (PFC) type model that expands the PFC methodol¬ 
ogy to a wider class of structurally complex crystal structures than previously possible. Specifically, our new 
approach allows for stabilization of graphene, as well as its coexistence with a disordered phase. It also pre¬ 
serves the ability to model the usual triangular and square lattices previously reported in 2D PFC studies. Our 
approach is guided by the formalism of classical field theory, wherein the the free energy functional is expanded 
to third order in PFC density correlations. It differs from previous PFC approaches in two main features. First, 
it utilizes a hard-sphere repulsion to describe two-point correlations. Second, and more important, is that it 
uses a rotationally invariant three-point correlation function that provides a unified way to control the formation 
of crystalline structures that can be described by a specific bond angle, such as graphene, triangular or square 
symmetries. Our new approach retains much of the computational simplicity of previous PFC models and al¬ 
lows for efficient simulation of nucleation and growth of polycrystalline 2D materials. In preparation for future 
applications, this paper details the mathematical derivation of the model and its equilibrium properties, and uses 
dynamical simulations to demonstrate defect structures produced by the model. 


I. INTRODUCTION 

Graphene is one of the most exciting new 2D material dis¬ 
covered. It has been found to exhibit interesting electrical 
Elll and mechanical properties mill. Crystals of graphene 
can be obtained by exfoliating cleaved graphite samples onto 
an oxidized silicon wafer to produce flakes of graphene 0. 
A more scalable method of obtaining graphene is through 
cvD gmi, a method that produces a polycrystalline material. 
While theoretically graphene can be about a hundred times 
stronger than steel, the properties of graphene realized in ex¬ 
periments typically reveal a wide variability, up to an order of 
magnitude from their theoretical predictions Eia. 

Variability in the strength properties of graphene, partic¬ 
ularly how these are related to the defect structure, still re¬ 
mains largely unexplored, particularly theoretically. Recent 
work suggests that it is linked to the defect microstructure at 
grain boundaries lIS UTOl . The topological defects in graphene 
typically take the form of periodic patterns of heptagonal and 
pentagonal disclinations, the patterning of which are dictated 
by the tessellation requirements of atoms in adjacent grains 
miHia. The complexity of forming and measuring graphene, 
however, make it challenging to experimentally isolate and ex¬ 
amine the role of specific defects and grain boundaries on the 
growth and properties of this material. 

Computational modelling can serve as a route for theo¬ 
retically understanding the difficult to measure properties of 
graphene. First principles studies are useful in examining 
the adsorption process of carbon onto metal surfaces dur¬ 
ing graphene formation m. Molecular dynamics (MD) 
studies of graphene have been successful at predicting the 
anisotropy of graphene morphologies on metal surfaces ca 
or the anergy of specific defect structures lfT3]l . On the con¬ 
tinuum scale, phase field models have been used to study how 
anisotropic diffusion of carbon on a surface can yield the for¬ 
mation of the dendritic graphene structures Qa. To date, 


there has not been a model that can address both the atomi¬ 
cally varying defect microstructures of graphene alongside its 
nucleation and diffusional growth kinetics from a disordered 
state on a surface. 

The phase field crystal (PFC) modelling approach is a 
promising approach for modelling many microstructure phe¬ 
nomena. The approach describes the thermodynamics and dy¬ 
namics of phase transformations through an atomically vary¬ 
ing order parameter field that is loosely connected to the 
atomic density field. Like traditional phase field (PF) mod¬ 
els, PFC models naturally capture most of the salient physics 
of nucleation, polycrystalline solidification, grain boundaries 
ifTTl - En and multi-component, multi-phase solidification ll2^ 
l25l . Unlike traditional PF models, PFC models also cap¬ 
ture, in the context of a single order parameter, elasticity and 
plasticity phenomena relevant to solid state processes such as 
dislocation source creation, dislocation stability 1^ l27l and 
creep ll^ . The most important feature of PFC-type models is 
that they incorporate the above phenomena from atomic to mi¬ 
cron length scales and over diffusional times scales, where the 
emergent properties of non-equilibrium phase transformations 
are typically manifested. PFC modelling has been used to elu¬ 
cidate phenomena ranging from grain boundary pre-melting 
mi im |30l to the nucleation pathways of defect-mediated 
nucleation of precipitates, the latter of which led to TEM 
and atomic probe experiments to validate the PFC predictions 

El- 

The original PFC model was predominately used for the 
study 2D triangular and 3D BCC crystal symmetries ifTTlI^ . 
Later models introduced multi-peaked two-point correlation 
kernels in the non-local part of the free energy that allowed 
for a simple yet robust manner to simulate most of the com¬ 
mon metallic crystal structures (2D Square, BCC,FCC,HCP) 
in phase transformations ll^ [34l . These so-called structural 
PFC (XPFC) models were later generalized to binary and 
multi-component (and multi-phase) alloys ll23l[T5ll . More re- 
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cently, a new multi-peaked two-point correlation was intro¬ 
duced to stabilize graphene and Kagome lattices, and a mor¬ 
phological phase diagram distinguishing the stability ranges 
of the two solid phases was explored numerically ll^ . 

This paper introduces a new structural PFC theory that 
breaks with the tradition of previous PFC models and ex¬ 
pands the free energy up to three-point correlations in the 
PFC density field. Unlike previous PFC theories, the two- 
point excess term is based on hard sphere-like interactions. It 
is shown that this allows for stabilization of triangular sym¬ 
metry in two dimensions. In this formalism, more complex 
crystal structures that are describable by a particular bond an¬ 
gle are stabilized using a new, rotationally invariant, three- 
point correlation function we introduce for the excess free 
energy. This term allows for the stabilization of triangular, 
square and graphene lattices in two dimensions. It is notewor¬ 
thy that beyond stabilizing the aforementioned structures, this 
formalism also allows for stable coexistence of these struc¬ 
tures with a disordered phase, a feature crucial for modelling 
nucleation and growth of polycrystalline 2D materials from a 
vapour or a disordered arrangement of atoms on a surface. In 
preparation for future applications, this paper highlights the 
derivation of our PFC free energy, examines its equilibrium 
properties and use dynamics to demonstrate defect structures 
produced in polycrystalline samples. 


II. NEW STRUCTURAL PFC MODEL 

The derivation of our model begins by defining the spatial 
PFC density field, p, of a species of atoms. From this, a di¬ 
mensionless density field is defined as n = {p — p)/p, where 
p is the reference density of a disordered phase around which 
a functional expansion of the free energy is carried out. Treat¬ 
ing the field n as an order parameter with which to describe 
microstructure variations, we expand the free energy of a crys¬ 
tallizing system as 

AF 

, rj.- = Fid[n] + Fex, 2 [n] + Fex,3[n] (1) 

kb-I p 

The term Fij_ is the ideal free energy, which ignores interac¬ 
tions. Its form here is given by 


energy and describes interactions at the level of three-point 
correlations in the density n. Its form, assuming translational 
invariance is given by 

Fex,3 = -^ Jn{r)Jcsir-r',r-r'')n{r')n{r'')dr'dr''dr 

(4) 

where C 3 is the three-point correlation function. The remain¬ 
der of this section will discuss the explicit forms of C 2 and 
C 3 , and their mathematical properties. 

A. Two-point correlations 

We define C 2 in Ff,x ,2 using a simple repulsive term. The 
length scale defined by this term will set the crystal lattice 
spacing. We define the circ function as a circular step function 
of unit radius and height centred on the origin, 

= { 0 : r > 1 

In two dimensions we then define C 2 as 

C'2(r) =--^circ ^ (6) 

’^’’0 VoJ 

where tq sets the cutoff for the repulsive term and R sets the 
magnitude of the repulsion. The normalization factor has been 
set such that in reciprocal space ( 72 ( 0 ) = —R, where (72 (k) 
denotes the Fourier transform of (72 (r). The form of (72 (r) is 
depicted schematically in Fig. [T] 

C 2 


0 
R 

FIG. 1: Two-point correlation function in real space. 
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( 2 ) 


It is convenient to simulate PFC models numerically in 
Fourier space. The Fourier transform of Eq. becomes 


where rj and x are dimensionless parameters to adjust the form 
of the ideal free energy. This form is a Landau expansion of 
the true ideal free energy. In what follows p = X = 1- The 
term 2 is the first term in the expansion of excess free 
energy, which incorporates two-point interactions. Its form is 
written as 

Fex ,2 = -^ jn{r)jc 2 {r-Y')n{r')dY'dr (3) 

where C 2 is the two-point correlation function. The term 
Fex, 3 [n] is the second term in the expansion of the excess free 


(72 (k) = (7) 

rok 

where Jm are the Bessel functions of the first kind. Figure 
shows a plot of (72 (k). Using Eq. 0, the excess free energy 
due to two-point correlations in Eq. Q can be written with the 
use of the convolution theorem as 

Fex,2 = -\ yn(r)F"^{(72(k)n(k)}dr (8) 

While there is no attraction between atoms at the level of 
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C2(k) 



FIG. 2: Plot of —2RJi{rok) /(rofc) in units of fc/ro. 


Lattice 

KiQq 

ro/ao 

Triangular 

4:TTjV3 

0.707854... 

Square 

27T 

0.81736... 

Graphene 

47 r/3 

1.22604... 


TABLE I: The ratio of ro to ao for various two dimensional crystal 
lattices. Ki is the reciprocal lattice vector of the first mode of a 
crystal structure. 


two-point correlations, the system can still undergo a phase 
transition and solidify when the density is great enough. The 
lattice constant oq is related in a nontrivial way to the cut¬ 
off length rg. Assuming for simplicity that most of the en¬ 
ergy is carried in the first reciprocal space mode (Ki) of the 
density, the energy of the system will be minimized when the 
first mode lies on the peak of —2RJi{rQk)/{rok), i.e., when 
Ki « 5.13562/ro. Table [I] gives the ratio tq/uq for various 
two dimensional lattices. 

Since the two-point correlation function has only one equi¬ 
librium distance and is completely isotropic it strongly favours 
a lattice with the highest packing fraction. In two dimen¬ 
sions this is the triangle phase, a result consistent with the 
classical result of hard sphere theory iJTll . To stabilize more 
structurally complex solid phases we must include either ad¬ 
ditional distances (as in XPFC models) or introduce a term 
which breaks the isotropy in interactions. Since we require 
our free energy to be rotationally invariant, this breaking of 
isotropy can only be relative to some local density configura¬ 
tion. This is not possible to achieve with a two-point corre¬ 
lation function - we must proceed to higher order and con¬ 
sider three-point correlation functions. This will permit us 
to favour particular relative angles between nearest neighbour 
atoms in order to produce crystal structures of interest, partic¬ 
ularly those of non-metals. 


B. Three-point correlations 


system). There is no such reduction in complexity available 
for Eq. Q to our knowledge. By comparison, it requires a 
prohibitive 0{N^) number of calculations, making it imprac¬ 
tical for most purposes. 


To remedy this problem we propose to separate C' 3 (r — 
r', r — r ) in the following manner 

C' 3 (r-r',r-r") = ^ - r')C7W(r - r") (9) 

i 

where the will be defined below. While this separation 
limits the possible forms of (73, it is sufficiently flexible to 
produce a wide variety of crystal structures, including the ones 
previously modelled by XPFC and similar 2D models. When 
Eq. (|^ is inserted into Eq. Q we obtain 

'W(r — r')n(r')dr'^ dr (10) 



Details on how to approach this term computationally will be 
discussed below. Here it should be apparent that in reducing 
the integration over r' and r" to one just over r' considerably 
reduces the computational complexity of the three-point term. 


We next define the Cl functions. We work in polar coor¬ 
dinates, and separate r into r and B according to 

Cf>{r,B) = Cr(r)c'^s\B) = (7^(r) cos(m6») (11) 

= Cr{r)cf'{B) = Cr{r) sm{m9) (12) 

Crir) =-^6{r - ao) (13) 

ZTTQjQ 

Where X is a parameter defining the strength of the interac¬ 
tion, og corresponds to the lattice spacing, and m defines bond 
order (discussed below) of the crystal phase. 


It may appear as though Eqs. -( [T^ break the isotropy 
of the free energy, but in fact they do not. While each term 
in the sum over i in Eq. ( [T0| ) exhibits angular dependence, the 
total sum is rotationally invariant. To see this, consider the 
sum in Eq. ( [T0| l and expand the square. This gives 

— r')C^'‘\r — r")n{r')n{r")dr'dr"'^ (14) 


Three-point correlations in the model are described by the 
excess energy in Eq. Q. This term is computationally expen¬ 
sive compared to the two-point excess term of Eq. Q. By use 
of the convolution theorem the two-point correlation can be 
computed by transforming to reciprocal space and multiply¬ 
ing point-wise. The computational complexity is therefore on 
the order of 0{N log N), the complexity of the fast Fourier 
transform (where N is the total number of grid points in the 


Defining ri = r — r', r 2 = r — r" and pulling the sum inside 
the integral gives 

= (i' 2 )^ n{r - ri)n(r - r 2 ) dri dr 2 

(15) 

Considering the sum in the brackets, changing to polar coor- 
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dinates ri —>■ (r^, 9i) gives 


III. EQUILIBRIUM PROPERTIES OE MODEL 


^CW(ri)C«(r2) =^C'«(ri,0i)C«(r2,02) (16) 

i i 

= Cr{ri)Cr{r2)Y.cf\B{)C^\e2) (17) 


(i) 

Inserting the Cg from Eqs. (Ill and (12 1 yields 


Y^Cf\v^)C^\v 2 ) = 

i 

Cr{ri)Cr{r 2 )^ cos(m0i) cos(m 02 ) 

+ sin(m0i) sin(TO6*2)| (18) 

Finally, noting the identity 

cos(0i) cos(02) + sin(0i) sin(02) = cos(02 — ^i) (19) 

gives 


^C«(ri)C«(r2) = a(ri)a(r2)cos(m(02-0i)) (20) 

i 

Since Eq. ( [20| ) depends only on the difference 62 —9i, the free 
energy remains isotropic. Moreover, by selecting the appro¬ 
priate values for m we can favour certain crystal bond angles 
over others. For example, m = 6 favours six-fold triangular 
crystals, to = 4 favours four-fold square crystals and to = 3 
favours three-fold graphene crystals. 


It will be useful below to have the Fourier transforms of 
(r, 9). Starting with (r, 9), we rewrite it as a multi¬ 
pole expansion. 


Ci^\r,9) 


X 

27rao 


S{r - ao) 



( 21 ) 


which transforms as 


(7(1)(It 9k) = 

' ( 22 ) 

pimOk I p-im9k 

= Xi^ -^-J™(fcao) (23) 

= Xi"^ cos{m9k)Jmikao) (24) 


where {k, 9k) are the polar coordinates in Fourier space. In 
Eq. (23 I we have used the fact that J_^{r) = (—1)"* Jm(r’). 
Proceeding similarly for 9) gives 


CP{k,9k) = Xi"^ sm{m9k)Jm{kao) (25) 


This section derives the equilibrium properties of our new 
structural PFC model defined by Eqs. 0-0 and with 
Eqs. ([TT|-([T3 ]i. In particular, we construct phase diagrams de¬ 
scribing solid-disorder coexistence for the case of three crystal 
systems; graphene, triangular and square. 

To describe the equilibrium properties of the model, we can 
expand the density of a crystal structure in a Fourier series; 

”■(!■)= (26) 

q 

where q are the reciprocal lattice vectors of the crystal struc¬ 
ture. The vectors q can be grouped according to their magni¬ 
tude; these groups are collectively referred to as modes, which 
are indexed by the integer k. The vectors within a mode are 
indexed by j, and the notation qkj is used to refer to a vector 
j of mode k. To simplify the equilibrium analysis, we assume 
that all amplitudes of a given mode are equal and real, i.e., 
(/fqj, J = (/)fc V j (see Figure|^. With these changes our expan- 



FIG. 3: Organization of amplitudes according to modes. Dots repre¬ 
sent reciprocal lattice vectors, circles connect vectors of equal mag¬ 
nitude, i.e. in the same mode. 

sion becomes 


(27) 

where we have truncated the expansion to N modes. 

To determine equilibrium properties of the model we use 
an approach well documented in numerous other PFC papers 
papers ll22]l . We begin by inserting Eq. ( |27| ) into the ideal and 
excess free energies and integrating over a unit cell. Since (pk 
are constants, the integration removes the dependency on r, 
leaving a free energy dependent on the amplitudes pk and the 
model parameters (R and X). The constant po corresponds 
to the average system density, while (pk, k > 0 , comprise a 
set of order parameters for the crystal phase. Minimizing the 
free energy numerically with respect to the (pk determines the 
phase and free energy of the system for a given {po, R, X}. 
When pk = 0 the system is in the disordered state. For non¬ 
zero pk the system is in a crystal state. Mapping out the con- 


k=o \ j 
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vex hull of the free energy as a function of ipo gives a common 
tangent line that defines coexistence values of (j>o correspond¬ 
ing to the solid and disordered phases. 


A. Triangular-disorder coexistence 

Using the above procedure with a five mode expansion 
for a 2D triangular phase, we produce a phase diagram for 
triangular-disorder coexistence. As noted above, the rejection 
term R is sufficient to produce a triangular phase without the 
need for three-point correlations. Figure shows the phase 
diagram in {t/fQ, R} space obtained by setting 2f = 0. 



B. Square-disorder coexistence 

To simulate square phases we require both two and three- 
point correlation functions in the free energy. To calculate 
the phase diagram for square-disorder coexistence, we use a 
five mode density expansion to describe the square phase, and 
set TO = 4. Figure shows a square-disorder phase diagram 
constructed with R = 5. 



FIG. 6: Square-disorder phase diagram using two and three-point 
correlations with m = 4 and R = 5. Here ro/ao = 0.81736. 


FIG. 4: Triangular-disorder phase diagram using only two-point 
correlations (X = 0). 


With R fixed we can also produce a triangular-disorder 
phase diagram in {(j)Q,X} space. Setting the three-point cor¬ 
relation function set to produce six-fold symmetry (to = 6), 
Figure l^shows such a phase diagram with R = 6. 



C. Graphene-disorder coexistence 


The graphene crystal structure can be described using a tri¬ 
angular lattice with a two atom basis. The primitive vectors of 
the triangular lattice are 


ao 


= VSx; 


V3 3 


while the basis atoms are located at 


. 73 1 

do = 2^’ 


^ 73 

di = --x+2y 


(28) 


(29) 


We must distinguish between graphene and triangular phases 
on the basis of the amplitudes (j)k of the density expansion. 
Assume we have a density expansion of a triangular phase tir 


as in Eq (27 1 . From this we can construct a graphene density 
field by 


FIG. 5: Triangular-disorder phase diagram using two and three- 
point correlations, with m = 6 and R — 6. Here ro/ao — 0.70785. 


ncir) = ^nH(r - d*) 


(30) 
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This gives 

; (31) 

(32) 

i 

where Sk,j are the structure factors for graphene. The basis 
vectors in Eq. ( |29l l have been selected such that the structure 
factors are independent of j, and so for our purposes = 
Sk- Computing the first four structure factors we find 

^0 = 2; 5i = -l; ^2 = 2; ^3 =-1; (33) 

Note that some modes have negative structure factors. This 
is in contrast to the structure factors of the triangular crystal 
structure, which are all unity. We can thus distinguish between 
the triangular and graphene phases by the signs of when 
(/)i < 0 we have a graphene phase and where > 0 we have 
a triangular phase. 

The presence of these negative amplitudes for graphene 
make it difficult to stabilize graphene-like structures using a 
two-point correlation function alone. This can be seen by in¬ 
serting the amplitude expansion for a triangular lattice into 
the two-point correlation term in Eq. ([^ and integrating over 
a unit cell, giving 


N 




fc =0 



FIG. 7: Graphene-disorder coexistence phase diagram, with m = 3 
and R = 6. Here ro/ag = 1.2259. 


feet structures at grain boundaries. 

The density n represents a conserved order parameter. As 
a result, its evolution is described by model B dynamics for 
conserved fields ll^ . This gives 


dn 

dt 


= m„v2 



(35) 


Fex,2 — — 2 (^2(9o)</’o + 6(72 (gi);/)! + 6(72 (<72)</>2 + ■ • ■ 


(34) 


Since all amplitudes in Eq. (34 1 are squared, the two-point cor¬ 
relation function alone cannot break the symmetry between 
positive and negative amplitudes, the latter of which are re¬ 
quired to stabilize graphene. Moreover, since the ideal free en¬ 
ergy is minimized by positive amplitudes, the triangular phase 
results. The three-point correlation function, however, results 
in terms of odd power which make it possible to minimize 
the free energy with negative amplitudes in order to produce a 
graphene phase. 

Using a five mode triangular density expansion in the free 
energy with to = 3 and applying the common tangent con¬ 
struction yields a graphene-disorder coexistence phase dia¬ 
gram in {00:7f} space. Eigure|^ shows such a phase diagram 
for the case R = 6. As expected, it was found that 0i < 0 for 
the ordered region of the phase diagram. 


IV. DYNAMICAL SIMULATIONS 


Where M is an effective mobility that sets the scale of the 
diffusional dynamics of n. Equation (35 i should also have a 
noise source to subsume the role of thermal fluctuations. In 
this work noise will be ignored. 


To simulate Eq. ( (T5| l, it is instructive to compile the various 
terms in the variational 6F/6n. The variational of the ideal 
term is straight forward. The variational of the two-point ex¬ 
cess term in Eq. Q becomes 


^Fex,2 

dn 



r')n(r')dr' = —C 2 * n 


(36) 


Where * indicates convolution. Terms such as this may be 
efficiently computed in reciprocal space by use of the convo¬ 
lution theorem. 

The three-point excess term of Eq. ( [T0) i is more complex. 
Eor each term in the summation we expand the square. 



r'')n{r)n{r')n{r'') dr' dr" dr 


(37) 


In this section we demonstrate the dynamical stability of 
the phase coexistence predicted by the phase diagrams of 
the previous section. We also demonstrate the robustness of 
our model to simulate nucleation, growth and formation or 
polycrystalline graphene, square and triangular crystal phases 
from a disordered phase. Eor the case of graphene, we also 
demonstrate the formation of experimentally relevant of de¬ 


We find the functional derivative using the well known for¬ 
mula 

f SF 

F[n + Sn] — F[n] = Sn^ dr (38) 

J on 

Applying this to Eq. ([T0|i and discarding terms of order 
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0(^{6n)^) we are left with three terms. 


Fex,3bT- + H - Fex,3M = 

— ^ Jc^'^\r — r')C^'‘\r — r")Sn{r)n{r')n{r")dr'dr"dr 

— - yC'^*^(r — r')C'^*^(r — r")n{r)Sn{r')n{r") dr' dr" dr 

— ^ Jc^''\r — r')C^^\r — r")n{r)n{r')Sn{r")dr'dr"dr 

(39) 


Swapping r and r' in the second term and r and r" in the third 
gives; 

FexA''"'+ - ^ex,3bA = 

~\ j (r - r")n(r')n(r") 
+ C^''\r' — r)C'f*^(r' — r")n(r')n(r") 

+ C^'\r" — r')C\''\r" — r)n(r')n(r")^ dr' dr" dr 

(40) 


Combining the last two terms (by swapping r' and r" in the 
third term) gives 

FexA''"'+ - Fex,3M = 

j6n{r)j (cA (r - r')CA (r - r")n(r')n(r") 

+ 2CA{F — ~ r")n(r')n(r")^ dr' dr"dr 

(41) 


Finally it is noted from Eqs. (Ill and (12 1 that r) = 

(-l)™(Ci*^(r). This leads to 


5F, 


\i) 


ere,3 


6n 


* n? 


-2{-ircfK[n-{cA*n)]) 

(42) 


The first term can be computed by performing the convolution 
in reciprocal space and then returning to real space in order to 
compute the square. The second term can be computed in 
several steps; performing the inner convolution in reciprocal 
space; returning to real space for the multiplication by n; and 
hnally transforming to reciprocal space once again to compute 
the outer convolution. 


A. Polycrystalline 2D materials, defects and coexistence 

We simulate the growth of triangular, square and graphene 
phases by choosing parameters corresponding to the solid re¬ 
gion of the phase diagram and initializing the system with 
gaussian noise. The system subsequently solidifies into a 
polycrystalline solid. Figure shows early and late time 



FIG. 8: Density fields of triangular (a and b), square (c and d) and 
graphene (e and/) growth, showing early (left) and late (right) times 
during solidification. Systems are initialized with gaussian density 
fluctuations. For all three systems /o = 0.3 and values of ro/oo 
match those in Table|^ For triangular, R = 7 and X = 0. For square, 
R — Q and X~^ = 0.5. For graphene, R = 6 and X~^ — 0.4. 


frames for the density of triangular, square and graphene sys¬ 
tems under crystallization. 

The defect structures which emerge along the graphene 
grain boundaries are noteworthy. Closer inspection of grain 
boundaries such as those in Figure]^) reveals that where mis¬ 
aligned grains impinge the grain boundary is resolved into a 
line of so-called 5-7 (pentagons and heptagons) defects. Fig¬ 
ure highlights these 5-7 defects. These defects are in ex¬ 
cellent agreement with the structures seen experimentally in 
polycrystalline graphene HOl . Like these experimentally de¬ 
termined structures, our simulated grain boundary consists of 
an aperiodic line of 5-7 defects. 

To demonstrate dynamic coexistence between the graphene 
and disordered phases, a 2000 x 100 system with periodic 
boundary conditions was seeded with a large initial slab of 
graphene and allowed to reach equilibrium with the disor¬ 
dered phase (Figure [T^g)). We set X~^ = 0.5 and R = 6. 
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FIG. 9: Comparison of simulated and experimentally determined 
defect structures of polycrystalline graphene. The defect structure of 
Figure[^) is highlighted in (a). The grain boundary is resolved by 
a line of 5-7 defect structures. These defect structures match those 
found experimentally in polycrystalline graphene membranes grown 
by chemical vapour deposition (CVD) l40l . (b) shows an atomic res¬ 
olution transmission electron microscope (TEM) image of one such 
graphene membrane; the defect structure is highlighted in (c). (b) 
and (c) reprinted by permission from Macmillan Publishers Ltd: Na¬ 
ture ESI, copyright 2011. 


Our approach differers from previous ones in two major 
ways. First of all, it treats two-point correlations more for¬ 
mally through the use of hard-sphere interactions. As a re¬ 
sult, the crystallography of structurally more complex phases 
than 2D triangular must be described in a unified way through 
the new rotationally invariant 3-point correlation introduced 
in this work. We showed that the form of our three-point cor¬ 
relation is rotationally invariant and robust enough to capture 
all crystal structures described through a single bond angle. 

After deriving the mathematical details of our new model, 
we calculated its equilibrium properties. We then used dy¬ 
namical simulations to illustrate the growth of poly crystalline 
graphene and other solids, and dynamical coexistence of 
graphene with a disordered phase. We also compared the 
defect structures generated at grain boundaries against corre¬ 
sponding results from the experimental literature, finding ex¬ 
cellent agreement. 



This approach, where the slab extends through the system 
traverse to the long dimension, negates the effects of curva¬ 
ture on the equilibrium coexistence densities since the order- 
disorder interface is a straight line. Figure [T0|fc) shows the 
smoothed density across the interface in the longitudinal di¬ 
rection. It can be seen that the equilibrium coexistence densi¬ 
ties in the ordered and disordered regions closely match those 
of the phase diagram in Figure]^ 


V. CONCLUSIONS 

This paper introduced the formalism of a new structural 
PFC theory that is truncated at three-point density correla¬ 
tions in the excess free energy. This approach makes it possi¬ 
ble to simulate microstructural evolution in metallic and non- 
metallic materials, as well as their coexistence with a dis¬ 
ordered phase. Among the most important novel materials 
that can be studied with our new formalism is polycrystalline 
graphene. 


FIG. 10: Simulation of coexistence between the ordered and disor¬ 
dered phases of graphene. Density field n(r) of the equilibrium inter¬ 
face between phases shown in (a). Smoothed average density along 
the longitudinal axis depicted in (b). X~^ = 0.5, R — 6. Average 
densities of 0.057 and 0.134 in the disordered and ordered phases 
respectively match closely the theoretical values from the phase dia¬ 
gram in Figure]^ 

It is expected that the structural PFC formalism introduced 
here is easily amenable to the recent formalism of Ref. Il25l . 
whereby the addition of an additional long-wavelength in¬ 
teraction energy term can be used to bring the pressure (P), 
X and 00 axes simultaneously under control. Similarly, our 
model is extendable to multiple components. These additions 
and their subsequent applications will be presented in future 
papers. 
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